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SUBROUTINE  READ 

DIMENSION  X(1000)»  YllOOO  t  BETAIlOUO),  ZETA'iiOOOl*  itlOOOJ 

COMMON  N.CCOE.  IL.  IU.A,EPS.GAMMA,>;,Y,?,E  lA.ZEl  A,S»BE  TAO.XI  .STOP 

l.MXSTP 

BMACHF(A)  s  SQRTFI A»A-1. > 

1  FORMAT(6I3.4£13.0J 

2  FORMAT{4Elb.0/£15.B) 

3  FORMATIEIO.OJ  ^  . 

4  FORMAT ( 43H1ERR0R  STOP.  INITIAL  FLOW  IS  SUBSONIC.  M  «  F12.6) 

READ  INPUT  TAPE  5 , 1 ,M. I  CODE . I U . I L . MAX .MXSTP , EPS .GAMMA , FREEM. ANGLE 
READ  INPUT  TAPE  S. 3. STOP 
EPS  =  EPS/IjO. 

N  =  M-1 

CODE  -  I  CODE 

A  =  MAX  O  D  W 

IF  (GAMMA)  lu.  10. 15  — 

10  GAMMA  =1.4  ""  "  M  I 

15  IF  (10-4)31.20.31  III  v.9?»bD4  |1 

20  IF  (FREEM-1.)  1000.  30.  30  i  .r  4  t.  j  i 

DLO  WRITE  OUTPUT  TAPE  6.  4.  FREEM  llnTTTCnn-all  \J  15111/ 

99  CALL  EXIT  “  B 

GO  TO  99  DDC-IKft  o 

3u  BETAO  *  BMACHF ( FREEM) 

XI  =  TANF (ANGLE) 

GO  TO  32 

31  BETA0=0. 

X1»0. 

32  DO  100  I»liN 

100  READ  INPUT  TAPE  5 . 2  .X  (  I  )  »  Y  (  I  )  t  t>E TA  (  1  )  .2t  T A  (  I  >  .S (  1  ) 

N=M-f  1 

READ  INPUT  TAPE  5 .2  »X  ( N)  . Y  I  f.)  1  BETA (  N )  »4t T A ( N  )  »S( N ) 

RETURN 

END 
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1,  Introduction 

In  an  article  published  in  the  Journal  of  the  Society  for 
Industrial  and  Applied  Mathematics  (Ref  Cl]),  the  writer  formulated  the 
method  of  characteristics  for  the  supersonic  isoenergetic  flow  of  an 
idesd  gas  specifically  for  coding  on  a  high  speed  digital  computer.  By 
the  introduction  of  variables  3  =  -  1  and  ^  =  tan  ©,  where  M 

is  the  local  Mach  number  and  0  is  the  local  flow  angle,  the  coefficients 
of  the  characteristic  equations  were  simplified  to  algebraic  rational 
functions,  thereby  appreciably  reducing  the  required  computing  time.  This 
document  describes  the  program  coded  for  the  IBM  7090  electronic  computer 
for  the  method  of  Ref  [1],  By  means  of  this  program,  it  is  possible  to 
compute  the  stxialTy-symmetric  or  two  dimensions^.  Mach  net  for  supersonic 
flows  over  pointed  bodies  with  attached  shocks,  through  channels  consisting 
of  two  rigid  walls,  and  through  nozzles  with  an  axis  of  symmetry.  The 
computation  of  flows  with  free  streamline  boundaries  is  also  possible,  if 
no  shocks  occur  within  the  flow, 

2,  Description  of  the  general  procedure  for  computing  the  Mach  net. 

The  complete  progrsun  uses  suiy  initisd.  data  line  which  does  not 
coincide  with  a  Mach  line  to  compute  the  location  and  flow  parameters  for 
the  next  set  of  points  downstream  of  this  initial  line.  To  start  the 
computations,  the  flow  data  along  this  initial  line  must  be  provided  and 
must  include  the  data  for  two  boundary  points  at  the  ends  of  the  line. 
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The  acheme  of  storing  the  points  Is  Illustrated  In  the  drawing  of  Fig.  1. 
Note  that  the  lower  boundary  point  data  Is  denoted  by  the  subscript  1. 
Interior  points  of  the  Mach  net  are  computed  by  a  subroutine  denoted 
by  ITER  (with  ITER  1).  By  means  of  the  CHAR  (k)  subroutine,  ITER  for 
odd  numbered  data  lines  takes  the  top  boundary  point  denoted  by  N  +  1 
and  the  next  lower  point,  N  -  1,  computes  the  next  downstream  point, 
and  stores  It  In  the  location  for  the  point  N.  For  other  points  on  the 
Initial  line,  ITER  takes  the  consecutive  points  j  and  j  +  1  ,  computes 
a  third  point,  and  stores  the  data  in  the  location  for  the  point  j  +  1. 
This  procedure  is  continued  down  the  line  of  data.  Note  that  the  Initial 
line  data  is  cancelled  out  except  for  the  boundary  points  and  therefore 
each  set  of  computatiors  must  be  stored  on  the  output  tape  for  later  off 
line  printing.  Since  the  upper  boundary  point  may  be  needed  in  the 
calculation  of  the  next  upper  boundary  point,  the  Initial  line  is  loaded 
to  skip  the  locations  for  the  point  N,  Thus,  the  second  line  of  data 
(and  all  even  numbered  lines  thereafter)  contains  data  for  N  -  1  points 
located  for  points  2  consecutively  to  N,  All  odd  numbered  data  lines 
start  with  the  point  1,  run  consecutively  to  N  -  1,  but  skip  N  with 
the  upper  boundary  point  in  the  location  for  N  +  1, 

When  even  numbered  data  lines  are  used  to  compute  the  next  downstream 
data  lines,  the  upper  and  lower  boundary  points  are  computed  first  and 
are  stored  in  locations  for  the  points  N  +  1  and  1,  respectively.  All 
remaining  points  are  computed  by  using  data  points  J  and  j  +  1  with  the 
new  data  stored  in  the  point  j  locations. 
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3.  Designation  of  Types  of  Flow. 

The  program  computes  the  Mach  net  for  either  isentropio  or 
non-isentropic  flow  and  for  either  axially-symraetric  or  two  dimensional 
flow  according  to  a  specified  code.  This  flow  code  is  an  integer  from 
1  to  4  denoted  in  the  progrsun  by  ICODE  or  CODE  depending  upon  the  form 
used  by  the  program.  The  interpretation  of  these  codes  is  given  in 
Table  I. 

TABLE  I 

ICODE  TYPE  OF  FLOW 

1  Axially  symmetric  with  variable  entropy 

2  Axially  symmetric  and  isentropic 

3  Two  dimensional  with  variable  entropy 

4  Two  dimensional  and  isentropic 

4.  Designation  of  upper  and  lower  boundary  points. 

The  selections  of  the  subroutines  for  the  upper  and  lower  boundaries 
are  made  by  code  numbers  denoted  by  ID  and  IL,  respectively.  The 
types  of  lower  boundaries  provided  for  are  the  axis  of  symmetry,  fixed 
wall,  and  the  free  streamline  bordering  still  air.  For  the  upper  boundary, 
the  fix  wsdl,  the  free  streaunline,  and  the  inclined  shock  with  a  uniform 
flow  ahead  of  it  may  be  chosen.  The  seune  code  numbers  are  used  for  both 
the  upper  and  lower  boundary  and  are  given  in  Table  2. 
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TABLE  2 


IL  or  m 
1 
2 

3 

4 


BOUNDAST 
Axis  of  symmetry 
Fixed  wall 
Free  streamline 

Shock  preceded  by  a  uniform  stream. 


Note  that  codes  IL  =  4  and  lU  =  1  are  Inadmissible  codes.  When 
these  codes  are  sensed  by  the  progreun,  the  calculations  are  stopped, 
a  statement  of  the  error  is  written  on  the  output  tape,  and 
calculations  are  terminated. 


5.  Preparation  of  input  data  cards. 

The  first  input  data  card  contains  the  code  numbers,  the  number  of 
points  and  other  parameters  associated  with  the  flow  field.  The  format 
of  this  card  is  expressed  in  Fortran  language  by 


FORMAT  (613,  4E13.0) 


The  order  of  the  data  on  the  first  card  is 


1.  M. 

2.  ICODE. 

3.  IP. 

4.  IL, 

5.  MAX. 


The  number  of  points  on  the  initial  data  line. 

A  number  from  1  to  4  designating  the  type  of  flow  and 
it  is  determined  according  to  Table  1  in  Section  3> 
Upper  boundary  code  according  to  Table  2. 

Lower  boundary  code  according  to  Table  2. 

Maximum  number  of  iterations  to  be  calculated  when  the 
iteration  does  not  converge  to  the  accuracy  specified. 


i 
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6. 


7. 

8. 

9. 

DO. 


MXSTP.  The  fflaxlmum  number  of  steps  (or  data  lines)  to  be 

computed.  This  must  be  an  even  number,  If  the  last 
data  line  Is  to  be  punched  on  cards  as  the  Initial 
line  for  continuing  the  Mach  net  at  a  later  run  time, 
EPS,  Maximum  percent  difference  between  the  values  of  each 
quantity  from  consecutive  Iterations, 

GAMMA.  The  value  of  the  ratio  of  specific  heats  for  the  gas. 
If  space  is  left  blank,  GAMMA  Is  set  equal  to  1.4. 
FRBEM.  The  free  stream  Mach  number  ahead  of  the  shock  wave. 
(Not  required  unless  lU  =  4). 

ANGLE  The  angle  of  inclination  of  the  shock  at  the  initial 
line  upper  boundary  point  for  code  ID  =  4, 


The  second  data  card  with  the  format  (lElO.O)  contains  STOP,  the 
maximum  value  of  X  to  which  the  Mach  net  is  to  be  computed.  If 
X  =  STOP  is  not  attained  in  the  number  of  steps  specified  by  MXSTP,  the 
program  writes  this  information  on  the  tape  to  be  printed  at  the  end  of  the 
calculations  and  punches  the  last  set  of  data  on  cards  for  continuing  the 
Mach  net  farther  downstream  in  a  later  run. 


The  cards  for  the  initial  line  are  loaded  in  the  machine  consecutively 
beginning  with  the  lower  boundary  point  data.  The  data  for  each  point  is 
punched  on  two  cards  according  to  the  format  (4E15.8/E15.8) .  The  order 
of  the  data  is 

1.  X  Axial  coordinate  of  the  point. 

2.  T  Vertical  coordinate  of  the  point. 

3.  BETA  The  loced  cotangent  of  Mach  angle,  p  =  -  1,  where  M 


is  the  Mach  number, 


7 


4.  ZSTA  The  local  tangent  of  the  flow  angle. 

5.  S  The  local  value  of  the  entropy.  (S  s  0  If  the 

flow  is  Isentroplc) . 

When  the  boundary  code  numbers  IL  =  2,  or  lU  =  2,  or  both  are 
prescribed,  then  FORTRAN  statements  must  be  Included  for  the  y 
coordinates  and  the  slope  of  the  boundary  curves  as  functions  of  the 
axial  variable  x  . 

The  subroutines  for  the  coordinates  and  slope,^  WALL(X,K)  and 
WALPR(X,K)  are  listed  on  page  56  of  Appendix  III  for  the  upper  wall 
given  by  the  parabola 


y  =  1  >  xVs 

and  for  the  lower  wall,  given  by 

y  =  0.1  +  0.05004l66x  -  0.010010i<2x^ 

The  BOUND  (M)  subroutine  supplies  the  code  K  =  0  for  computing  the 
upper  wall  and  K  =  1  for  the  lower  wall. 
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APPENDIX  I t  Interaction  of  a  Mach  Line  vdth  a  Shock  Wave. 


Since  the  procedure  for  computing  the  interaction  of  a  Mach  line 
with  a  shock  wave  described  in  Ref.[l]  does  not  always  converge  we  have 
chosen  a  different  method  of  iteration,  which  we  shall  explain  in  seme 
detail.  Let  the  quantity  at  the  previous  shock  point  be  denoted  by 
the  subscript  2  and  the  nearest  interior  point  by  1  (see  Fig.  2).  The 
first  estimate  for  the  location  of  the  new  point  on  the  shock  is  given 

by 

=  (^2  ~  *  ^*1  " 

y3  =  Yl  +  A^(x^  -  x^)  (2) 


where  is  the  slope  of  the  shock  at  the  point  (2)  and  is  defined 
in  Ref,[l],  The  flow  variables  at  the  new  point  must  satisfy  the 
compatability  relation 


(3) 

+  H^^y^  “  y^^Xi/yj^  +  i/y^/2.  =  o. 

As  shown  in  section  9  of  Ref  [1],  P3*  functions  of  an 

the  Mach  number  ahead  of  the  shock.  Thus  the  left  hand  side  of  Eq.(3) 

may  be  regarded  as  a  function  of  ®®y  then  that  value  of 

which  satisfies  Eq,(3)  can  be  determined  by  Newton's  method,  i.e., 
by  iterating 


=  53  -  f(j:3)/f'(i:3) 


! 


(4) 
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To 


evaluate  f(llj)  and  f'(JIj)  we  use  the  formulas  of 


Sj  from  Ref.[l].  These  are  written  in  the  following  form: 

=  gC-  1  +  0^(5^  +  Ds^jOo  ^  +  2)(c^i:|  -  1)] 

C  d  -  1  (C  d  +  2) 

B  =  log[-^-2 5  ]  +  Y  log[ - 5 — -  ] 


g(l  + 


(y  +  i)i:|(3o  +  1) 


where 


g  =  (y  +  1)/(y  -  1), 


=  (y  -  l)pQ  +  Y  +  1, 

C2  =  +  2(Pq  +  1), 


(5) 

(6) 

(7) 


(8) 


Cj  =  C^/(y  -  1)  +  gpQ  . 


Differentiating  Eqs.  (5)  and  (6)  logarithmically  and  also  Eq,  (7),  we  obtain 


for  the  derivatives  with  respect  to  5^, 


1  +  C, 


pi  =  (p^  +  — I - 2  2 

3  3  Y  - 1  3  +  2)  (1  +  d)(cj:‘^  - 1) 


VP  (9) 


’3"'l'’5  ■  '  -3''  3^ 

=  2(2p2  -  +  c^)  -  C5  /  ?5  (10> 

1  +  C  - 

s'  =  2i:j - - 2 - 5 - T—^ -  3  (11) 


^  ^  *  4^  4^^l4  * 


With  Eqs.  (9),  (10),  and  (11),  is  given  by 


f(j;j)  =  E^C^  -  F^p^  -  Q^s^ 


(12) 
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In  the  process  of  iterating  Eq,  (4),  estimates  of  S^,  and 
are  obtained  as  well  as  of  5^.  The  iteration  is  continued  until  the 
value  of  converges.  The  sequence  of  calculations  is  repeated  by 

replacing  the  coefficients  of  Eq,(3)  and  (12)  by  their  average  values 
between  the  points  1  and  3  and  between  the  points  2  and  3.  For 
example  is  replaced  by  (F^  +  Fj)/2  in  Eqs.  (3)  and  (12)  and 

and  are  replaced  by  and  +  A^)/2  in  Eqs.d)  and 

(2).  An  alternate  procedure  is  to  evaluate  A^  ,  E(C),  F(p),  etc,,  for 
the  quantities  ( Pj^  +  P^)/2  and  (^^^  +  t^)/2. 

The  program  as  listed  on  page  37  Appendix  III  has  an  additional 
feature  which  evens  out  the  M'ich  net  in  the  region  of  the  shock.  When  the 
Mach  net  is  started  from  an  initial  line  including  a  point  on  the  shock, 
the  next  point  on  the  shock  wave  will  be  farther  downstream  than  the  other 
points  on  that  line  of  data  (compare  points  3  and  3'  of  Fig,  2).  This 
caul  be  remedied  by  computing  the  next  point  on  the  shock  from  a  new  interior 
point  found  by  interpolating  linearly  between  the  two  points  A  eind  B 
eilong  the  Mach  line  AB  in  Fig.  3.  Lines  BC  and  ED  are  chosen  to 
have  the  slope  of  the  Mach  line  at  B  given  by 


A 

n 


Pn^n  ^  ^ 
Pn-^n 


(13) 


The  equations  of  line  DE  and  AB  are 


y  -  A^(x  -  xp 


y  -  yn.l  =  ^l^’'  -  Vl^ 


(14) 


(15) 


u 


reapectlvely,  where 


D.  a  (y  ,  -  y  )/(x  ,  -  X  ), 
1  'n+l  *'n  n+1  n 


(16) 


Choosing  ~  ***  determine  y^  by  the 

intersection  of  line  DE  with  the  shock  line  given  by 


yn>i  “  yf  =  Vi " 


(17) 


where  J[  is  the  slope  of  the  shock  at  •  Solving 

Eqs.  (14)  and  (15)  for  x  and  y  and  substituting  the  values  of  x^ 


and  y^  yields 


*  =  Vl  *  2(?  -  A^)(x^  - 


y  =  yn+1  ^  -  Vi>- 


(18) 


By  the  point  (x,y),  the  segment  AB  is  divided  into  segments  whose 
lengths  have  the  ratio  D2  /  (  1  -  D^)  where 


°2  =  ^Vl  -  “  ^n>' 


Hence,  the  values  of  b  at  the  point  E  are  given  by 

equations  of  the  form 

f  =  -  ‘>2>- 

The  next  point  along  the  shock  is  computed  by  using  the  values  at  the 
interpolated  point  in  place  of  the  values  at  interior 

point  denoted  by  the  subscript  1  in  the  foregoing  discussion. 
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APPENDIX  II;  The  Design  of  Perfect  Nozzles 

Since  the  program  as  described  in  the  main  part  of  this  document 
is  capable  of  computing  a  symmetric  isentropic  flow  with  a  fixed  wall, 
it  is  possible  to  make  use  of  it  in  the  design  of  both  two  dimensional 
and  axially  symmetric  perfect  nozzles.  However,  an  additional  program 
is  required  to  compute  the  transition  region  which  produces  a  uniform 
supersonic  flow. 

A  method  for  designing  supersonic  axially  symmetric  nozzles  consisting 
of  a  circular  throat  region,  of  a  conical  expanding  section,  euid  of  a  final 
transition  to  uniform  flow  (see  Fig. 4),  is  described  in  Ref. [2].  The 
assumption  that  the  flow  field  in  the  coniceil  region  is  approximated  by  the 
supersonic  source  leads  to  a  set  of  equations  for  finding  the  transition 
boundary  Mach  line.  For  a  nozzle  with  a  cone  section  of  half  angle 
and  length  L  and  with  a  throat  radius  I^,  the  coordinates  of  the 
transition  Mach  line  are  given  parametrically  by 

x/h  =  (t  cos  ©  -  Oj^)/2TgBin(©^/2) 

+  (Hj/h)8in  ©^  +  (I/h)cos  ©  (19) 

y/h  =  (T/'''g)ein  ©/2sin(©^/2)  . 

The  parameter  ©  is  the  local  flow  angle,  h  is  the  throat  radius,  and 
T  is  related  to  the  local  value  of  P  by 

T  =  (1  +  p^/k)^V(P^  +  1)^^  (20) 


where 


k  =  (y  +  1)/(y  -  1) 
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The  subacripts  A  and  E  denote  quantities  evaluated  at  the  points 
A  and  E  of  Fig.  4,  respectively.  The  quantities  p  and  C  =  tan  0 
along  the  Mach  line  are  functionally  related  by 

6  =  (uj,  -  («))/2  (21) 

where  u  is  the  Prandtl-Meyer  angle, 

0)  =  tan~^(p/^)  -  tan~^P  (22) 

which  is  tabulated  in  Ref. [3]  page  504  and  following.  With  Eqs.  (19) 
through  (22),  the  values  of  x,y,p,  and  C  for  points  on  the  transition 
boundary  Mach  line  are  determined  parametrically  by  a  choice  of  P. 

Dividing  the  range  of  p^  to  Pg  into  N,  say,  intervale,  we  obtain  a 
smoothly  spaced  set  of  points  along  the  boundary  Mach  line. 

By  means  of  the  foregoing  procedure,  or  by  means  of  the  program 
described  in  the  main  section  of  the  document,  we  can  establish  values  of 
P,4,x,  and  y  at  N  +  1  regularly  spaced  points  along  the  Mach  line  AE 
in  Fig.  (4),  Since,  we  specify  uniform  flow  at  the  exit  end  of  the  nozzle, 
the  values  p  =  Pg  and  C  =  0  are  known  along  Mach  line  EC  and  x 
and  y  are  then  related  by 

V  =  ■  ^E  • 

With  conditions  specified  along  the  two  Mach  lines  AE  and  EC,  the  flow 
is  uniquely  determined  in  a  region  AECB  bounded  by  the  four  Mach  lines 
shown  in  Fig,  4. 
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By  working  upward  step  by  step  along  the  Mach  line  AE  from  a 
point  on  the  line  £C,  we  calculate  the  next  Mach  line  of  the  same 
family  (see  Fig. 5).  This  is  easily  accomplished  by  using,  for  example, 

the  simple  program  suggested  in  Table  3,  page  19,  together  with  the 
subroutine  ITER  and  with  the  STREM  routine  programmed  by 
Malcolm  Gray  [4]  and  listed  in  Table  4,  page  2o. 

The  subroutine  STREM  utilizes  a  method  of  polynomial  interpolation 
to  find  the  streamline  point  on  the  next  downstream  Mach  line.  Three 
points  are  chosen  so  that  the  y  coordinate  of  the  highest  point  is 
greater  than  the  value  of  y  at  the  last  streamline  point,  x^,  y^. 
Denoting  the  points  on  the  Mach  line  by  subscripts  1,2,  and  3,  (see  Pig, 6) 
and  applying  Lagranges'  interpolation  formula  yield  the  following  equations 
for  y(x),  C(x)»  and  P(x)  along  the  Mach  line} 

y(x)  =  yj^Lj^(x)  +  y2L2(x)  +  yjLj(x) 
i;(x)  =  Cj^Lj^(x)  +  +  CjLj(x) 

P(x)  =  Pj^I.^(x)  +  +  PjLj(x) 

where 

Lj^(x)  =  (x  -  X2)(x  -  Xj)/(xj^  -  x^Xxj^  -  Xj) 

L2(x)  =  (x  -  Xj^)(x  -  Xj)/(x2  -  "  *3^ 

L^(x)  =  (x  -  Xj^)(x  -  x^)/(.x^  -  x^)(xj  -  X2) 

The  line  with  the  slope  (^q  +  ^(x))/2,  representing  the  streamline 
sgment  is  given  by 

F(x)  =  y  -  y^  -  [C(x)  +  Cq](x  -  Xq)/2  =  0 
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The  solution  by  Newton's  method  of  this  last  equation  with  y  -  y(x) 
along  the  Mach  line  yields  the  next  downstreaun  point  on  the  boundary 
streamline  of  the  transition  region. 

The  program  is  capable  of  using  np  to  50  points  along  the  Mach  line. 
The  variables  X(I),  Y(I),  Z(l),  B(l)  denote  the  Cartesian-coordinates 

t 

of  the  points  on  the  Mach  line  and  their  corresponding  values  of  C 
and  p,  respectively.  The  quantities  XYBZO(l),  I  =  1,2,3»  and  4, 
denote  the  values  of  X,Y,P,  and  Ct  respectively,  for  the  last  stream¬ 
line  point,  and  XYBZl(I),  the  corresponding  values  for  the  point  to  be 
computed  on  the  streamline.  The  symbols  XY(1)  and  XY(3)  represent 
the  coordinates  x,y  of  the  point  E  in  Fig.  4,  and  XY(2)  and 
XY(4),  for  the  point  C. 

The  code  number  IDONE  =  1  when  the  program  is  entered.  At  exit, 
the  code  number  may  be  changed  and  then  has  the  following  interpretations: 

IDONE  INTERPRETATION 

1  Last  streamline  y  coordinate  is  less  than 

the  Isurgest  vsdue  of  y  on  the  Mach  line. 

Next  point  of  the  boundary  is  successfully 
computed. 

2  Last  streamline  y  coordinate  is  greater  than 

the  largest  value  of  y  on  the  Mach  line.  New 
point  is  not  determined. 
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IDONE  IWrERPRETATICW 

3  Smallest  y  on  Mach  line  is  greater  than 

y  of  previous  streamline  point.  This 
indicates  that  the  next  streamline  point  is 
on  Mach  line  EC  of  Fig,  4  and  the  point  is 
determined  under  this  assumption. 

4  The  iteration  for  x  did  not  converge  to  the 

required  relative  accuracy  EPSIL  in  10 
iterations. 

The  number  N  is  set  equal  to  the  number  of  the  lowest  point  for  which 
the  y  coordinate  is  greater  than  the  previously  computed  boundary  point. 
The  transition  region  program  in  Table  3  uses  this  feature  of  STREM  to 
reduce  systematically  the  number  of  points  computed  for  each  successive 
Mach  line. 

To  compute  two  dimensional  nozzles,  relations  similar  to  Bqs,  (19) 
through  (22)  can  be  derived  for  the  boundary  Mach  line, or  with  a  known 
boundary  for  the  throat  region  of  the  nozzle  the  Mach  net  ceui  be 
calculated  by  the  program  described  in  Sections  1  to  5.  By  assuming  that 
the  lines  of  const£uit  Mach  number  are  circles  and  are  derived  from  a  source 
flow,  one  cam  obtain  the  initial  line  data  (see  Ref. [3],  and  Fig.  7).  The 
Mach  number  is  found  from  the  ratio  of  the  area  of  the  initial  line  surface 
to  the  throat  area.  Calculation  of  the  transition  region  for  a  two 
dimensional  nozzle  is  much  simpler  since  the  region  consists  of  one  family 
of  straight  Mach  lines  along  which  the  flow  properties  are  constant.  3y 
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approximating  the  streamline  boundary  between  two  Mach  lines  (see  Fig.  8) 
by  a  straight  line  whose  slope  is  the  average  of  the  slopes  on  the  two 
Mach  lines,  we  can  determine  the  location  of  next  boundary  point  by 
solving  two  simultaneous  linear  equations  (see  Ref, [53).  In  this  way, 
we  obtain 

*  ^  ^^n+l  “  ^0  ^^n  *  ^n+l^*©'^^  "  \+l\+l^/^^^n  *  ~  *n+l^ 


y  =  y 


n+1 


+  A  ,(x-x  ,) 
n+1  n+1 


where 


A  ,  =  O  ,  + 

n+1  ^n+l^n+1 


l)/(0 


n+1 


■*  ^n+1^ 


and  the  subscripts  denote  the  values  of  the  quantities  at  the  corresponding 
points  shown  in  Fig,  8.  A  program  suggested  for  computing  the  Prandtl- 
Heyer  region  is  given  in  Table  5.  Note  that  the  data  for  the  boundary 
Mach  line  are  read  beginning  at  the  upper  and  upstream  end  of  the  line.  The 
point  1  must  also  be  a  boundary  point  on  the  fixed  wall. 
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TABLE  3 


C  TRANSITION  REGION  AXIALLY  SYMMETRIC  NOZZLE 

DIMENSION  X(50) *Y(50) *BETA( 50 ) tZETAI 30 ) »XYBZ0( 4 ) iXYBZl ( 4 ) »XY ( 4 ) 

10  READ  INPUT  TAPE  5t20»  M*N*EPS*DX*GA 
20  FORMAT  (2I3*3E16.8) 

READ  INPUT  TAPE  5  *  30 1 C  XYBZO  (  H  *  1  >  1  *4  )  •  (  X  Y  (  I  )  •  1  i  ,4  >  ,  ( x  (  1  )  i  Y  (  I  )  tBET 
1A(  I  )>ZETA(  1  )  tl^^l.N) 

30  FORMAT  (4E16.8) 

XODE  =  2. 

IDONE«l 
S  =  0. 

XE*X( 1 ) 

NN»N-1 
DO  60  I-1*M 
X(1)=X(1)+DX 
Y(1)=(X«1)-XE)/BETA(1) 

DO  40  JsltNN 

40  CALL  IT£R(BETA( J+1 )  .BETA( J ) »ZETA( J+1 » .ZETA ( J ) .S »S »X ( J+1 ) tXC J» . 

1  Y( J+1  )  .Y ( J) .KODEtN»EPS.GA.BETA( J+1 ) .ZETA (J  +  1 ) . S ( J  +  1 ) .X ( J+1 ) . Y ( J  +  1 
1 )  ) 

N=NN+1 

CALL  STREM(N. IDONE.EPS.X.Y.BETA.ZETA.XYBZO.XYBZI .XY  > 

NN»N-1 

WRITE  OUTPUT  TAPE  6 . 30 » ( XYBZ 1 ( I ) . 1 = 1 .4 ) 

DO  50  K*1.4 
50  XY8Z0(KI=XYBZ1 (K) 

GO  TO  (60.70.10.90) .lOONE 
60  CONTINUE 
GO  TO  10 

70  WRITE  OUTPUT  TAPE  6.80 

80  FORMAT (42H  ERROR  RETURN  FROM  STREM — K  GREATER  THAN  Nl 
GO  TO  10 

90  WRITE  OUTPUT  TAPE  6.100 

100  FORMAT(52H  ERROR  RETURN  FROM  STREM-- 1 TERAT I  ON  DID  NOT  CONVERGE) 

GO  TO  10 
END 


INPDT  DATA 


N  =  Number  of  Intervals  along  straight  Mach  line 

BC  in  ng.  4. 

N  =  Number  of  points  along  transition  boundary 
Mach  line  A£  in  Fig.  4. 

DX  =  Interval  length  in  x  for  points  along  EC. 

QA  =  y,  the  ratio  of  specific  heats  for  the  gas. 
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TABLE  4 


SUBROUTINE  STREM  (N.  IDONE.  EPSIL»  X.  Yt  B>  Z.  XYBZO»  XYBZl.  XY) 
DIMENSION  X(50)t  Y(50)»  B(50).  2(50).  XYBZ0(4).  XYBZ1(4).  XY(4) 
INDIC  *  1 
K  »  1 

1  IF  (XYBZ0(2)  -  Y(K)  )  2.  4.  8 

2  IF  (K  -  1)  18.  18.  3 

3  <  =  K  +  1 
GO  TO  20 

4  IF  (XYBZ0(4))  5.  5.  7 

5  XYBZl(l)  =  X(K) 

XYBZ1(2)  =  Y(K) 

XYBZIO)  =  B(K) 

XYBZ1(4)  »  0. 

IF  (K  -  1)  6.  6.  21 

7  IF  (K  -  2)  18.  18.  8 

8  <  »  K  +  1 

IF  <K  -  N)  1.  9.  19 

9  IF  (XYBZ0(4)  -  Z(N))  20.  20.  19 

19  IDONE  =  2 
GO  TO  17 

20  N  =  K 

X1MX2  =  X(K)  -  X(K-l) 

X1MX3  =  X(K)  -  XiK-2) 

X2MX3  =  X(K-1 )  -  X(K-2) 

GIDEN  »  X1MX2»X1MX3 
G2DEN  =  -X1MX2*X2MX3 
G3DEN  =  X1MX3»X2MX3 
H  =(Y(K)  -  Y(K-l) )/XlMX2 

XIOLD  =  (XY8Z0(2)  -  Y(K)  +  H*X(IC.)  -  XYBZO  (  1  )  *X  YBZO  (  4  )  )  /  (  H  -  XYBZ0( 
14)  ) 

LOOP  =  1 

10  TERMl  =  XIOLD  -  X(K) 

TERM2  =  XIOLD  -  X(K-l) 

TERM3  »  XIOLD  -  X(K-2) 

G1  =  (TERM2»TERM3)/G1DEN 
G2  =  (TERM1*TERM3)/G2DEN 
G3  =  (TERM1#TERM2)/G3DEN 

XYBZ1(2)  =  Y(K)#G1  +  Y(K-1>»G2  +  Y(K-2)*G3 
XYBZ1«4)  =  Z(K)»G1  +  Z(K-l)*G2  +  Z(K-2)»G3 
GO  TO  (11.  16)  .  INDIC 

11  GIPRI  =  (TERM2  +  TERM3)/G1DEN 
G2PRI  =  (TERMl  +  TERM3)/G2DEN 
G3PRI  =  (TERMl  +  TERM2)/G3DEN 

YPRIM  =  Y(IC)»G1PRI  +  Y(K-1)*G2PRI  +  Y(K-2)*G3PRI 
H  =  XIOLD  -  XYBZ0( 1  ) 

ZPRIM  =  (Z(K)»G1PRI  +  2(K-1)»G2PRI  +  Z ( <-2 ) ♦G3PR I ) *H/2. 0 
ZFRAC  *  (XYBZ0(4)  +  XYBZ 1 ( 4 ) ) /2 .0 

XINEW  =  XIOLD  -  (XYBZ1(2)  -  XYBZ0(2)  -  H*ZFRAC ) / ( YPK I M  -  ZPRIM  - 
IZFRAC) 

IF  ((ABSF(X1NEW  -  XIOLD) /XINEW)  -  EPSIL)  12.  12.  13 
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12  INDIC  »  2 
XYB21(1)  «  XINEW 
GO  TO  15 

13  LOOP  »  LOOP  ♦  1 

IF  (LOOP  -  lOJ  15*  15*  14 

14  lOONE  >  4 
GO  TO  17 

15  XIOLD  «  XINEW 
GO  TO  10 

16  XYBZ1(3)  =  B(K)*G1  +  B(IC-1)*G2  +  B(K-2)»G3 

17  RETURN 

18  SLOPE  =  (XY(4)  -  XY(3) )/(XY(2)  -  XY(1)) 

XY5  =  XYBZ0(4)/2.0 

XYBZKl)  »  (XY(3)  -  SLOPE»XY(l)  -  XYBZ0(2)  +  XY5*XYbZ0(  H  ) 
1/(XY5  -  SLOPE) 

XYBZ1(2)  s  XY5*(XYBZ1(1)  -  XYBZO(l))  +  XYBZ0(2) 

XYBZ1(3)  s  B(l) 

XYBZ1(4)  »  0.0 
6  IDONE  =  3 
GO  TO  17 
21  N  *  K 
GO  TO  17 
END 


TABLE  5 


C  PRANDTL-MEYER  REGION 

DIMENSION  X(50) »Y(50) .BETA(50) .ZETA ( 50 ) » XB ( 50 ) » YB ( 50 ) 
APLMIF(AtC.M)  =  (A*C+(-l.)**M)/(A-(-l.)**M*C) 

10  READ  INPUT  TAPE  5*20. N 
20  FORMAT  (  13  ) 

READ  INPUT  TAPE  5 • 30* ( X ( I ) * Y ( I ) .BETA ( I > »ZETA ( I ) *  I>1*N) 

30  FORMAT  (4F12.8) 

XB(1)»  X(l) 

YB(1)«Y(1) 

DO  40  Is2*N 

A1  =  APLMIFI BETA< I ) *ZETA( I ) *0» 

Z1  =  (ZETA(I)  +  ZETA( I-l) )/2. 

XB(I)  =  (YBd)  -  YB(I-l)  +  Z1«XB(I)  -  A1*X  ( I  )  )  /  ( Zl-Al ) 

40  YBd)  »  Yd)  +  A1«(XB(I)  -  XCD) 

WRITE  OUTPUT  TAPE  6*30*  (XBd  )  *YBd  )  *BETA(  I )  *ZETAd  )  *  I«1*N) 

GO  TO  10 

END 


IHPOT  DATA 


N  =  Number  of  points  adong  boundary  Mach  line 
of  Prandtl-Meyer  Region. 
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APPENDIX  III  FQgraAM  Lj-«tlng  of  Sttbroutlnee 

for  the  ttethod  of  Characteristics 


C  MAIN  PROGRAM 


DIMENSION  X(IOOO).  YdOOO).  BETA(IOOO)*  ZETA(IOOO).  S(IOOO) 

COMMON  N.CODE»lL.IU.A*EPS*GAMMA»X»Y.BETA»ZETAfS.BETAOtXI .STOP 
1 .MXSTP 
80  CALL  READ 
N=N 

MXSTP=MXSTP 
CALL  WRITEl 
CALL  WRITE2(1) 

DO  100  1=1.  MXSTP 

IF  ( (-1 )**I )  10.  10.  20 
10  K  =  1 
GO  TO  30 
20  K  =  0 
30  CALL  CHAR(K) 

CALL  WRITE  2(1+1) 

DO  60  J=l.  N 

IF  (X(J)-STOP)  100.  60.  60 
60  CONTINUE 
GO  TO  80 
100  CONTINUE 

70  WRITE  OUTPUT  TAPE  6.  1 

1  FORMAT (64H  THE  LOWER  LIMIT  WAS  NOT  REACHED  IN  THE  MAXIMUM  NUMBER  0 
IF  STEPS.) 

M  =  N-2 

DO  150  1  =  1. M 

150  WRITE  OUTPUT  TAPE  14 .200.X ( I ) . Y ( I ) .BET A ( I ) .ZET A (  I  ) .S (  I  ) 

200  FORMAT (4E15.8/E15. 8  ) 

WRITE  OUTPUT  T APE14 .200 .X ( N ). Y ( N )» BETA ( N ) .ZET A( N ) .S ( N ) 

GO  TO  80 
END 
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SUBROUTINE  READ 


DIMENSION  XdOOOIt  YdOOOl*  BETAdOOOIt  ZETAdOOO)*  SdOOO) 
COMMON  NtCOOE*IL»IUtA*EPS*GAMMA*X*y*BETA*ZETAtS*BETAO*XI (STOP 
IfMXSTP 

BMACHF(A)  »  SP"TF(A»A-1.) 

1  FORMAT(6I3(4E13.0) 

3  FORMAT(EIO.O) 

READ  INPUT  TAPE  5 ( 1 (M t ICODE t lUt IL (MAX (MXSTP (EPS (GAMMAtFREEM t ANC 
READ  INPUT  TAPE  5(3(STOP 
EPS  *  EPS/100. 

N»M-1 

CODE  »  ICODE 
A  -  MAX 

IF  (GAMMA)10(10(15 
10  GAMMA  «  1.4 
15  IF  (IU-4)30(20(30 
20  IF  (FREEM-1.)  1000(  30(  30 
30  BETAO  »  BMACHF(FREEM) 

XI  s  TANFIANGLE) 

DO  100  I=1(N 

100  READ  INPUT  TAPE  5 (2 (X ( 1) . Y ( 1) ( BET A ( 1 ) (ZET A ( I ) (S d ) 

2  FORMAT(4E15.8/E15.8) 

N=M+1 

READ  INPUT  TAPE  5 ( 2 (X ( N ) t Y ( N ) (BET A ( N ) (ZETA ( N ) ( S ( N ) 

RETURN 

1000  WRITE  OUTPUT  TAPE  6(  4(  FREEM 
CALL  EXIT 

4  FORMAT (43H1ERR0R  STOP.  INITIAL  FLOW  IS  SUBSONIC.  M  »  F12.6) 

END 
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SUBROUTINE  WRITE  ) 


DIMENSION  X(1000)>  Y(IOOO)*  BETA(IOOO)*  ZETA(IOOO)*  S(IOOQ) 

COMMON  NfC0DE»IL»IU*AiEPS*6AMMA*X*Y»BETAiZETA*StB0*XI tSTOP* 

IMXSTP 
J  *  CODE 

WRITE  OUTPUT  TAPE  6.  1 
GO  TO  (I0t20»30.40) »  lU 

10  WRITE  OUTPUT  TAPE  6»  2 
GO  TO  50 

20  WRITE  OUTPUT  TAPE  6.  3 
GO  TO  50 

30  WRITE  OUTPUT  TAPE  6t  4 
GO  TO  50 

40  BMACH  s  SQRTF(B0*B0+1. ) 

WRITE  OUTPUT  TAPE  6.  5.  BMACH 
50  GO  TO  (60*70.80 ) t  IL 
60  WRITE  OUTPUT  TAPE  6.  6 
GO  TO  90 

70  WRITE  OUTPUT  TAPE  6.  7 
GO  TO  90 

80  WRITE  OUTPUT  TAPE  6.  8 
90  GO  TO  (100.100.110.110).  J 
100  WRITE  OUTPUT  TAPE  6.  9 
GO  TO  120 

110  WRITE  OUTPUT  TAPE  6.  11 
120  GO  TO  (130.140.130.140).  J 
130  WRITE  OUTPUT  TAPE  6.  12 
GO  TO  150 

140  WRITE  OUTPUT  TAPE  6.  13 
150  PCTER  =  EPS»100. 

WRITE  OUTPUT  TAPE  6. 14.PCTER.MXSTP.STOP 
RETURN 

1  F0RMAT(1H1///39X.42HCHARACTERISTIC  SOLUTION  OF  SUPERSONIC  FLOW) 

2  F0RMAT(///35X.27HILLEGAL  UPPER  BOUNDARY  CODE) 

3  FORMAT (///35X.27riUPPER  BOUNDARY  FIXED  WALL) 

4  FORMAT (///35X.28HUPPER  BOUNDARY  FREE  STREAM) 

5  FORMAT (///35X.42HUPPER  BOUNDARY  SHOCK  INITIAL  MACH  NO.  .E15.8) 

6  FORMAT (//35X.33HLOWER  BOUNDARY  AXIS  OF  SYMMETRY) 

7  FORMAT(//35X.27HLOWER  BOUNDARY  FIXED  WALL) 

8  FORMAT! //35X.28HLOWER  BOUNDARY  FREE  STREAM) 

9  FORMAT(//35X.22HAXIALLY  SYMMETRIC  FLOW) 

11  FORMAT(//35X.20HTWO  DIMENSIONAL  FLOW) 

12  FORMAT (//35X.12H I SOENERGETIC) 

13  FORMAT(//35X»10HISENTROPIC) 

14  FORMAT(//35X.12HNO  MORE  THAN.Fll. 6.33H  PERCENT  ERROR  WILL  BE  INTRO 
1DUCED.//35X.12HIN  ONE  STEP .// /35X . 34HTHE  PROCESS  WILL  BE  STOPPED  A 
2FTER  .I4.19H  STEPS  IF  THE  LOWER .// 35X . 1 IHL IMI T .  X  =  F12.6.17H.  IS 
3NOT  REACHED. ) 

END 


[ 
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SUBROUTINE  WRITE2(K) 


DIMENSION  X(IOOO)*  YIIOOO)*  BETA(IOOO)*  ZETA(lOOO)*  S(IOOO) 

COMMON  N»  CODE*  IL*  IU»  A*  EPS*  GAMMA*  X*  Y*  BETA*  ZETA*  S*  BO*  XI 

I  » 

10  WRITE  OUTPUT  TAPE  6*  1*  K 
M*N-2 

DO  1000  J»1*M 
II=I+J-1 
750  L=J 

lOUO  WRITE  OUTPUT  TAPE  6 *2 *L *X ( 1 1 ) » Y ( 1 1 ) *BET A ( 1 1  )  *ZETA ( 1 1)  *S (  1 1) 

IF  (<-l)»*K)  2000*2000*2050 
2000  L»N-1 

WRITE  OUTPUT  TAPE  6 *2 *L *X( N ) * Y ( N ) *  BETA ( N ) *ZETA ( N ) *S ( N ) 

IF  (IU-4)  2050*2020*2050 
2020  ANGLE  »  ATANFIXI  ) 

WRITE  OUTPUT  TAPE  6*  4*  ANGLE 
2050  RETURN 

1  FORMAT (1H1/39X.42HCHARACTERISTIC  SOLUTION  OF  SUPERSONIC  FL0W*///54 

1X*8HSTEP  N0.*I4////.16X,1HX*22X.1HY.18X.9HBETA  *13X*11HZETA 

2  *14X*7HENTR0PY*/120X) 

2  F0RMAT(120X/I5*5E23.8) 

4  FORMAT (120X/41H  THE  SHOCK  ANGLE  AT  THE  BOUNDARY  POINT  IS*E20.8) 

END 
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SUBROUTINE  CHAR(IC) 


DIMENSION  X(IOOO)*  Y(IOOO).  BETA(1000)»  ZETAdOOO)*  S(IOOO) 
COMMON  N»  CODE*  LOW.  HIGH.  A.  EPS.  GAMMA.  X.  Y.  BETA.  ZETA.  S 
M  »  N-2 

IF  (K)  1000.  10.  100 
10  CALL  BOUNDd  ) 

CALL  B0UND(2) 

IF  DIVIDE  CHECK  40.  50 
40  WRITE  OUTPUT  TAPE  6.  2 

2  FORMAT(//23H  DIVIDE  CHECK  IN  BOUND  I 
50  CONTINUE 

DO  15  J=2.  M 
15  CALL  ITERl ( J. ( J+1 I .J) 

IF  DIVIDE  CHECK  20.  30 
20  WRITE  OUTPUT  TAPE  6.  3 

3  FORMAT(23H  DIVIDE  CHECK  IN  ITER  ) 

30  CONTINUE 

RETURN 

100  CALL  ITER1(N-2.N.N-1) 

IF  DIVIDE  CHECK  110.  120 
110  WRITE  OUTPUT  TAPE  6.  3 
120  CONTINUE 

DO  150  J=2.  M 
I  =  M+2-J 

150  CALL  ITERKI-l.I.I) 

IF  DIVIDE  CHECK  160.170 
160  WRITE  OUTPUT  TAPE  6.  3 
170  CONTINUE 
RETURN 

1000  CALL  EXIT 
END 
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SUBROUTINE  ITERl  (IiJ*K) 


DIMENSION  X( 1000) »BETA( 1000) *ZETA( 1000 ) »S ( 1000 ) tY(lOOO) 

COMMON  N*COOE*IL*IU*A»EPS*GAMMA*X*Y»BETA»ZETA»S 

MAX«A 

)C0DE»C0DE 

CALL  I  TER (BETA! I ) iBETAI J) »ZETA( I ) *ZETA( J) tS( I)»S(J)*X(I)»X(J)*Y(I) 
1 tY( J) tKOOEtMAXiEPStGAMMAtBETAIK) tZETAOC) »S(K) »X(K) »Y(K) ) 

IFISENSE. LIGHT  4)10*30 
10  WRITE  OUTPUT  TAPE  6*20.I*J.K 

20  FORMATOOH  ITER  DOES  NOT  CONVERGE  FOR  1*15. 4H  J*I5.4H  K«I5) 

3C  CONTINUE 
RETURN 
END 
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INTERIOR  POINT  ITERATION  ROUTINE 


SUBROUTINE  I TER(B1 *B2*Z1»Z2 »S1 »S2 *X1 »X2 * Y1 (YZ * ICODEf MAX»EPS»GtBN3* 
lZN3*SN3iXNafYN3) 

FUNCTIONS 

EFNF(Z)  »  1.0/(1.0+Z»*2) 

APLF<B»Z)  ■  (B*Z+1.0)/ (B-Z» 

AMIF(B»Z)  =  (B*Z-1.0)/(B+Z) 

FFNF ( B ) » ( 2 . 0*B*»2 ) / ( ( 8**2+l .0 ) * ( ( G-1 .0  >  *B»*2+G+ 1 .0 ) I 
HYPF(B*Z)»Z/(B-Z) 

HYMF(B.Z)*Z/(B+Z) 

GFNF(B»»B/(G*(G-1.0>»(B«*2+1.0)  ) 

ROUTINE 

GO  TO  (20*20*10«10  )«ICOOE 
10  HYPl  =  0.0 
HYM2  *  0.0 
ICY  *  3 
GO  TO  50 

20  IF  (Yl)  40.  30.  40 
30  ICY  =  2 
GO  TO  50 
40  ICY  »  1 
50  IFT  =  1 
BOLD  >  0.0 
20LD  »  0.0 
SOLD*0.0 
XOLO  =  0.0 
YOLO  *  0.0 

GO  TO  (60.70.90).  ICY 
60  HYP1=HYPF(B1 .Z1 ) 

60  TO  80 

70  HYPl  =  Z2/(B1*Y2) 

80  HYM2=HYMF(B2.Z2) 

90  SHYPl  a  HYPl 
SHYM2  «  HYM2 
API  a  APLF(Bl.Zl) 

SAPl  a  API 

AM2  »  AMIF(B2.Z2) 

SAM2  a  aM2 
El  a  EFNF(ZI) 

SEl  *  El 
E2  a  EFNF(Z2) 

SE2  *  E2 
FlaFFNFIBl ) 

SFlaFl 

F2«FFNF(B2) 

SF2aF2 

GO  TO  (94.92.94.92) -(CODE 
92  G1»0.0 
G2a0.0 
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GO  TO  96 
94  G1=GFNF(B1) 

G2»GFNF(B2) 

96  SG1»G1 
SG2-G2 

Z3«(Z1+Z2  )*0.5 
S3«0.0 

DO  250  1=1.  MAX 

X3  =  (Y2-Y1+AP1*X1-AM2*X2J/(AP1-AM2) 

Y3  =  Y2+AM2*( X3-X2 » 

GO  TO  (98.100.98.100) .ICODE 

98  S3=S2+< (Y3-Y2+Z3*(X2-X3) )*(S1-S2) ) / ( Yl-Y2+Z3» ( X2-X1 ) ) 
100  GO  TO  (102.102.112.112 ) .ICODE 
102  GO  TO( 110.104.112 )  .ICY 
104  GO  TO  (106.108) .IFT 
106  IFT=2 

FJ3=(X3-X1)»Z2/(B1»Y2) 

GO  TO  114 

108  FJ3=(HYP3/Y3+Z2/(B1»Y2) )*(X3-X1)*0.6 
GO  TO  114 

110  FJ3=HYP1*(X3-X1  )/Yl 
GO  TO  114 
112  FJ3=0.0 

114  FJ=E1#Z1-F1*B1+G1*(S3-S1)-FJ3 

FK=E2*Z2+F2»B2-G2»(S3-S2)+HYM2»(X3-X2 )/Y2 
Z3=(FJ*F2+FK»F1 )/ (E1»F2+E2*F1 ) 

B3=(FI<;-E2*Z3  )  /F2 
IF(B3) 130.135.130 

130  IF  (ABSF( (B3-BOLD)/B3)-EPS)  135.  190.  190 
13b  IF(Z3)140.145.140 

14C  IF  (ABSF( (Z3-Z0LD)/Z3)-EPS)  145.  190.  190 
145  IF(S3)150.155.150 

150  IF  (ABSF( (S3-S0LD)/S3)-EPSI  155.  190.  190 
155  IF(X3)160.165.160 

160  IF  (ABSF( (X3-XOLO)/X3)-EPS)  165.  190.  190 
165  IF(Y3)170.180.170 

170  IF  (ABSF( (Y3-YOLD)/Y3)-EPS)  180.  190.  190 
180  BN3=B3 
ZN3=Z3 
SN3=S3 
XN3=X3 
YN3=Y3 
RETURN 

190  BOLD  =  B3 
ZOLD  =  Z3 
SOLD  =  S3 
XOLD  =  X3 
YOLD  =  Y3 

API  =  (SAP1+APLF(B3.Z3) )*0.5 
AM2  =  (SAM2+AMIF(B3.Z3) )*0.5 
E3»EFNF(Z3) 

E1=(SE1+E3)*0.b 

E2=(SE2+E3)»0.5 

F3»FFNF(B3) 

F1«(SF1+F3)*0.5 
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F2*'(SF2+F3)*0.5 
GO  TO  (200*210*200»210) tICOOE 
200  G3>GFNF(B3) 

Gl>(SGl-i-G3)*  0.5 
G2-  (SG2+G3)*0.5 

210  GO  TO  (220*220i250t230) tlCOOE 
220  HYP3*HYPF(B3.Z3> 

GO  TO  (240*  230i250)>ICY 
230  HYP1»(HYP3/Y3+SHYP1 )*0.5 
GO  TO  242 

240  HYP1»(HYP3*(Y1/Y3»+SHYP1  )*0.5 

242  HYM2«( ( 23»Y2 ) / ( ( B3+Z3 J *Y3 ) +SHYM2 > »0 .5 

230  CONTINUE 

SENSE  LIGHT  4 
GO  TO  180 
END 


MOTE:  This  subroutine  differs  slightly  from  the  calculation 
procedure  outlined  In  paragraph  2,  page  88  of  Ref.  [1], 
In  place  of  the  term 

jH*(p,C)dy/y 

In  Eq.  (11)  of  Ref.  [1]  we  have  programmed 

JH*0,0«lx/y, 

where  ^  is  now  defined  by 

H*(P,a  =  C/(P  ♦  c). 

This  eliminates  the  singularity  in  H*(p,^)  for  Mach 
lines  with  zero  slope. 


SUBROUTINE  BOUND  (M) 


DIMENSION  SPACE(2) 

COMMON  SPACE* ILrIU 
<  «  M-1 

GO  TO  (100*200) *  M 
100  GO  TO  (1000*120*130*140)*  lU 

120  CALL  FIXED(K) 

IF  DIVIDE  CHECK  121*  125 

121  WRITE  OUTPUT  TAPE  6*  4 

4  FORMAT (//23H  DIVIDE  CHECK  IN  FIXED  ) 

125  IF(SENSE  LIGHT  4)126*128 

126  WRITE  OUTPUT  TAPE  6*127 

127  FORMAT(52H  MAXIMUM  NO.  OF  ITERATIONS  EXCEEDED  IN  FIXED  ROUTINE) 

128  CONTINUE 
RETURN 

130  CALL  FREE(K) 

IF  DIVIDE  CHECK  131*132 

131  WRITE  OUTPUT  TAPE  6*2 

2  FORMAT (//23H  DIVIDE  CHECK  IN  FREE  ) 

132  IF(SENSE  LIGHT  4)133*135 

133  WRITE  OUTPUT  TAPE  6*134 

134  FORMAT(51H  MAXIMUM  NO.  OF  ITERATIONS  EXCEEDED  IN  FREE  ROUTINE) 

135  CONTINUE 
RETURN 

140  CALL  SHOCK 

IF  DIVIDE  CHECK  141*  142 

141  WRITE  OUTPUT  TAPE  6*  5 

5  FORMAT (//23H  DIVIDE  CHECK  IN  SHOCK  ) 

142  IF(SENSE  LIGHT  3)143*145 

143  WRITE  OUTPUT  TAPE  6*144 

144  FORMAT (62H  MAX.  NO.  OF  ITERATIONS  EXCEEDED  IN  SHOCK  ROUTINE  COMPUT 
lING  XI ) 

145  IF(SENSE  LIGHT  4)146*148 

146  WRITE  OUTPUT  TAPE  6*147 

147  FORMAT (52H  MAXIMUM  NO.  OF  ITERATIONS  EXCEEDED  IN  SHOCK  ROUTINE) 

148  CONTINUE 
RETURN 

200  GO  TO  (210*120.130*2000)*  IL 
210  CALL  AXIS 

IF  DIVIDE  CHECK  220*230 
220  WRITE  OUTPUT  TAPE  6*3 

3  FORMAT (//23H  DIVIDE  CHECK  IN  AXIS  ) 

230  IFISENSE  LIGHT  4)231*233 

231  WRITE  OUTPUT  TAPE  6*232 

232  FORMAT(51H  MAXIMUM  NO.  OF  ITERATIONS  EXCEEDED  IN  AXIS  ROUTINE) 

233  CONTINUE 
RETURN 

1000  WRITE  OUTPUT  TAPE  6*  1*IU 
CALL  EXIT 

2000  WRITE  OUTPUT  TAPE  6*  1*IL 
CALL  EXIT 

1  FORMAT (67H  PROGRAM  STOPPED  DURING  BOUNDARY  DUE  TO  AN  IMPROPER  BOUN 
IDARY  CODE*  II) 

END 


34 


SUBROUTINE  AXIS 


DIMENSION  X(1000)«  ¥(1000)>  BETA(1000)t  ZETA(IOOO)*  S(IOOO) 

COMMON  N*  CODE*  LOw.  HIGH*  CLIM*  EPS*  GAMMA*  X*  Y*  BETA*  ZETA*  S 
APLMIF(A*B*N)  =  ( A*B+ ( -1 . ) »*N ) / ( A-B» ( -1 . ) »»N ) 

EIEZFIA)  =  l./(A**2+l.) 

FIFZFIA)  =  2.»A*»2/( (A»*2+l.»»< IGAMMA-1. )*A**2+GAMMA+1. » ) 

G1G2F(A)  =  A/ (GAMMA»(GAMMA-1. )#(A»*2+1. ) ) 

HPLMIF(A*B*N)  =  B/ ( A*B+ ( -1 . ) **N ) 

IF  (CODE-3.)  12*  11*  10 

10  G2  s  0. 

H2  =  0. 

GO  TO  25 

11  H2  »  0. 

G2  »  G1G2F(BETA(2 ) ) 

GO  TO  25 

12  IF  (COOE-1.)  1000*  20*  15 
15  G2  *  0. 

GO  TO  17 

20  G2  *  G1G2F{BETA(2 ) ) 

17  H2  =  HPLMlF(BETA(2)*ZtTA(2)*l) 

25  M  =  CLIM 

E2={E1E2F(ZETA(2) )+l. )/2. 

F2  =  F1F2F(BETA(2 ) ) 

DO  50  1=1*  M 

bNEW=(E2»Z£TA(2)-G2*(S(l)-S(2) ) - ( ri2-ZE T A ( 2 ) )/2.  )/F2  +  bETA(2  > 

IF  (ABSF( (BNEW-BETA( 1 )  )/BETA( 1 ) )-EPS)  ICO*  100*  30 
30  F2=F1F2F( (BNEW+BETA(2) )/2. ) 

BETAd)  =  BNEW 
IF  (G2)  35*  50*  35 
35  G2=G1G2F( (BNEW+BETA(2) )/2. ) 

5C  CONTINUE 

SENSE  LIGHT  4 

lUO  X(l)  =  X(2  )-2.*Y(2 )/( APLMIF(BETA( 2 ) .ZETA(2 ) *1 )-l./BNEW ) 

BETA( 1 )  =  BNEW 
RETURN 

1000  WRITE  OUTPUT  TAPE  6.  1*  CODE 

1  FORMAT (49H  PROGRAM  STOP  DUE  TO  AN  IMPROPER  FLOW  CODE  NUMBER*F8.2) 
CALL  EXIT 
END 


35 


SUBROUTINE  FIXED  (K) 


DIMENSION  X(1000)»  Y(1000)»  BETAIlOOOIt  ZETAHOOOIt  S(IOOO) 
COMMON  Nt  CODE*  LOMi  HIGH*  A*  EPS*  GAMMA*  X*  Y*  BETA*  ZETA*  S 
APLMIF(A*B*M>  »  ( A»B+ ( -1 . ) •*M ) / ( A- ( • » •♦M^B ) 

E1E2F(A)  =  l./(A**2+l.) 

F1F2F(A)  s  2. •A*»2/ (  (A»*2+l.)*( (GAMMA-1. )*A**2+GAMMA+1. ) > 
G1G2F(A)  s  A/(6AMMA*(GAMMA-1.)»(A»*2+1.I ) 

HPLMIF(A*B*M1  =  B/ ( A*B+(-l. )**MI 


IF 

(K)  1000*  10*  20 

10 

Jl 

=  N-1 

J3 

=  N 

GO 

TO  30 

20 

Jl 

=  2 

J3 

=  1 

30 

A1 

=  APLMIF(BETA( Jl) *ZETA(J1) *K) 

El 

=  E1E2F(ZETA(J1) ) 

FI 

=  F1F2F(BETA( Jl ) ) 

IF 

(CODE-2.)  80*  40*  50 

40 

Gl 

=  0. 

GO 

TO  85 

50 

IF 

(CODE-3.)  1000*  70*  60 

60 

Gl 

=  0. 

65 

HI 

=  0. 

GO 

TO  90 

70 

Gl 

=  G1G2F(BETA( Jl) ) 

GO 

TO  65 

80 

Gl 

=  G1G2F(BETA( Jl ) ) 

85 

HI 

=  HPLMIF(BETA(J1) *ZETA(J1) *K> 

90 

X(J3)  =  l.ElO 

Y(J3)  =  I.EIO 
BETA(J3)  =  l.ElO 
ZETA(J3)  =  l.ElO 
XNEW=X( Jl) 

LIM  =  A 

DO  200  1=1.  LIM 

XNE»i=(  Y(  Jl)-wALL(XNEw»<)+XNEW»WALPR(XNEW.lL)-X(  J1)*A1  )  /  (  WALPR (  XNEW * 
IK)  -All 

YNEw  =  WALL(XNEW.K) 

ZNEW=WALPR(XNEw,K) 

BNEW  =  (  (-1. )*»K»E1»(ZNEW-ZETA( Jl ) )-Gl*(St J3>-S( Jl ) )+Hl»( YNEW- 
lY(Jl) )/Y(Jl) )/Fl+BETA(Jl) 

IF  (ABSF( (X( J3)-XNEW)/XMEW)-EPS)  100*  100*  150 
100  IF  (ABSF( (Y( J3)-YNEW»/YNEW)-EPS)  110*  110*  150 
no  IF  (ABSr ( (ZETA( J3)-ZNEw)/ZNEW)-EP5)  120*  120*  150 
120  IF  (ABSF( (BETA( J3)-8NEW)/BNEW)-EPS)  210*  210*  150 
150  X(J3)  =  XNEW 
Y(J3)  =  YNEW 
BETA(J3)  =  BNEW 
ZETA(J3)  =  ZNEW 
B1=(BNEW  +  BETA( Jll )/2. 

Z1=(ZNEW  +  ZETA(Jl) )/2. 

A1=APLMIF(B1*Z1*K) 

FlsFlF2F(Bl ) 
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E1«E1E2F(Z1  I 
IF  (Gl)  160*  170»  160 
160  G1=G1G2F(B1  ) 

170  IF  <Hl)  180.  200.  180 

180  H1  =  {HPLMIF(B1.Z1.KJ  )»(2.*Y(J1) ) / ( Y ( J1 ) +YNEW  ) 

200  CONTINUE 

SENSE  LIGHT  4 
21C  RETURN 

1000  WRITE  OUTPUT  TAPE  6.  1.  CODE 
CALL  EXIT 

1  FORMAT (42rilERROR  EXIT  FROM  FIXED  SUBROUTINE.  CODE  =  F7.4) 
END 


FUNCTION  WALL  (X.K) 

IF  (K»  1000.  10.  20 
10  WALL  =  1.  +  X»X/5. 

RETURN 

20  WALL=  .1  .05004166*X  -  .01001042*X*X 

RETURN 

1000  WRITE  OUTPUT  TAPE  6.  1.  K 
CALL  EXIT 

1  FORMAT! 38H1ERROR  EXIT  FROM  WALL  SUBROUTINE.  K  =  .16) 
END 


FUNCTION  WALPR(X.K) 

IF  (1C)  1000.10.20 
10  WALPR  =  ,4»X 
RETURN 

20  WALPR  »  .05004166  -  .02002084*X 
RETURN 

1000  WRITE  OUTPUT  TAPE  6.1.)C 
CALL  EXIT 

1  FORMAT! 38H1ERROR  EXIT  FROM  WALL  SUBROUTINE.  K  =  .16) 
END 
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SUBROUTINE  SHOCK 


DIMENSION  X(1000)»  Y(IOOO).  BETAIIOOO)*  ZETA(1000)»  S(IOOO) 

COMMON  N»  CODE*  XYZ*  ZYX*  C*  EPS*  UAMMA*  X*  Y*  BETA*  ZtTA*  S*  bO* 
IXI 

APLMIF(A*B)  »  (A»B+1. )/(A-B) 

E1E2F(A)  =  l./(A*A+l.l 

F1F2F(A)  =  2. *A*A/( (A»A+1. )*I (GAMMA-1. )*A*A+(GAMMA+1. )) ) 

G1G2F(A)  =  A/(GAMMA#(GAMMA-1. )*(A*A+1. ) ) 

HPLMIF(A*B)  =  B/<A*B+1.> 

Cl  =  (GAMMA-1. )*B0*BC+GAMMA+1. 

C2  =  C1+2.»(B0*B0+1.) 

GA 1  =  ( G AMMA+ 1  .  )  /  (.GAMMA- 1 .  ) 

C3=  Cl/ (GAMMA  -  1.1  +  GAl»BO»BO 

X3  =  l.E+10 

Y3  =  l.E+10 

Z3  =  l.E+10 

B3  =  l.E+10 

S3  =  l.E+10 

SB=BETA(N-1 ) 

SZ=ZETA(N-1 ) 

SX=X(N-1 ) 

SY=Y(N-1 ) 

SXI=XI 

A1  =  APLMIF(BETA(N-1) *ZETA(N-1I ) 

Dl=  (Y(N-l)  -  Y(N) )/(X(N-l)  -  X(N)) 

XC=  X(N)+2.0»(Al-xn»(X(N-l  )-X(N)  )  /(Al-Dl) 

YC=Y(N-1 )+01»(XC-X(N-l) ) 

02=(XC-X(N) )/(X(N-l)-X(N) ) 

BETA(N-1 )=BETA(N-1 ) +D2+BETA ( N ) * ( 1 . -D2 ) 

ZETA(N-1 )=ZETA(N-1 )»D2+ZETA(N)«( 1.-D2) 

S(N-1)=S(N^1 )»D2+S(N)*(1.-D2) 

X(N-1 )=XC 
Y(N-1 )=YC 

El  =  E1E2F(ZETA(N-1 )  I 
FI  =  F1F2F(BETA(N-1 I ) 

A1  =  APLMIF(BETA(N-1» *ZETA(N-1 )  ) 

G1  =  G1G2F(BETA(N-1 ) ) 

IF  (COOE-3. )55. 35*35 
35  HI  =  0. 

GO  TO  60 

55  HI  =  HPLMIF( BETA(N-1 ) *ZETA(N-1 ) ) 

60  L=C 

DO  300  1=1*  L 

XNEW  =  (Y(N-1)-Y(N)+XI»X(N)-A1»X(N-1) )/(Xl-Al) 

YNEW  =  Y(N-1  )+Al»(XNEW-X(N-l) ) 

XI2  =  XI 
DO  70  J=l*  L 
XIS  =  XI2»XI2 
FDl  =  C1«XIS  +  2. 

FD2  =  XIS  +  1. 

FD3  =  C3*XIS  -  1.  , 

FDA  =  FDl  +  C2  -  2. 

SFl  *  2./(XIS»FDl» 
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SF2  =  (C3  +  l.)/(FD2»FD3) 

2NEW  =  2.*(B0»B0*XIS  -  1 . ) / ( X  1 2*FD4 ) 

SNEW  =  GAMMA#L0GF(FD1/ ( (GAMMA  +  1 . ) *X 1 S» ( BO»BO  +  !.))» 

1  +L0GF(FD3/(GA1*FD2) ) 

BNSU  =  GA1*(-1.  +  (C1*FD2*GA1»XIS»(B0»B0  +1 . ) ) / ( FD1*F03 ) ) 
BNEW  =  SaRTF(BNSQ) 

ZP  =  2.*(2.*B0*B0-C1*XI2*ZNEW)'/FL>4-ZNEW/XI2 
5P  =  2.*XI2*(SF2  -  GAMMA*SFl) 

BP  =  XI2«(BNSQ  +  GA1)*(SF1  -  SF2)/BNEw 

FXI  =  E1»(ZNEW  -  ZETA(N-1)1  -  F1»(BNEW  -  BETA(N-l)) 

1  -G1*(SNEW  -  S(N-l))  +  H1*2.*(YNEW  -  Y ( N-1 ) ) / ( YNEW+Y ( N- 1 ) ) 
FXIP  =  E1*ZP  -  Fl»dP  -  G1*SP 
XI3  =  XI2  -  FXI/FXIP 

IF  (ABSF( (XI3-XI2 )/XI3)-EPS>  100»100.65 
66  XI2  =  XI3 
70  CONTINUE 

SENSE  LIGHT  3 

100  IF  (ABSF( (X3-XNEw)/XNEW)-EPS)  110.  110.  200 
no  IF  (ABSF( (Y3-YNEW)/YNEw J-EPS)  120.  120.  200 
120  IF  (ABSF(  (B3-BNEW)/aNEKv)-EPS)  13C.  130.  200 
130  IF  (ABSF(  (Z3-ZNEw)/BNEW)-EPS)  1-40.  140.  200 
140  IF(SNEW) 146.150.146 

145  IF(ABSF(  (S3-SNEH)/SN£vk)-EPS>160.160.200 
150  X(N)  =  XNEW 
Y(N)  =  YNEW 
BETA(N)  =  BNEW 
ZETAIN)  =  ZNEW 
S(N)  =  SNEW 
XI  =  XI3 
X(N-1)=SX 
Y(N-1)=SY 
ZETA<N-1)=SZ 
BETA(N-1  )=SB 
RETURN 

200  X3  =  XNEW 
Y3  =  YNEW 
Z3  =  ZNEW 
B3  =  BNEW 
S3  =  SNEW 
XI=(SXI+X13)/2. 

B1=(BNEW+BETA(N-1 ) )/2. 

Z1=(ZNEW+ZETA(N-1 ) )/2. 

A1=APLMIF(B1  .Zl  ) 

E1=E1E2F(Z1  ) 

F1  =  F1F2F(B1  ) 

G1  =  G1&2F(81  ) 

IF  (HI)  216.  300.  216 
216  Hl=  HPLMlF(Bl.Zl) 

300  CONTINUE 

•SENSE  LIGHT  4 
GO  TO  150 
END 
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SUBROUTINE  FREE  (M) 


DIMENSION  XdOUOIi  TdOOOi*  BETAdOUOIf  ZETAdOOO)*  5(1000) 
COMMON  N>  COOEi  LOw«  HIGH*  A*  EPS*  GAMMA*  X*  Y*  BETA*  ZETA*  S 
APLMIF(A*B*X )  =  (A*B+(-l,)»»K)/(A-B*(-l.)**<) 

E1E2F(A)  =  l./(A«»2+l.) 

F1F2F(A)  =  2. *A»A/( 1A*A+1. )*( (GAMMA-1. )*A*A+(GAMMA+1,) ) ) 
G1G2F(A)  =  A/ (GAMMA*(GAMMA-1. )»(A»A+1. ) ) 

HPLMIF(A*B*K)  =  B/ ( A*B+ { -1 . ) ) 

IF  (M)  1000*  10*  20 
10  J1  =  N-1 
J2  =  N 
K  =  0 
GO  TO  30 
20  J1  =  2 
J2  =  1 
X  =  1 

30  B1=(BETA< J1 )+BtTA( J2) )/2. 

A1=APLMIF(B1*ZETA(J1) *X) 

El  =  E1E2F(2ETA( J1 ) ) 

F1=F1F2F(61 ) 

IF  (CODE-2.)  80*  AO.  60 
40  Gl  =  0. 

GO  TO  90 

50  IF  (CODE-3.)  1000.  60*  70 
61/  )I1  =  0. 

Gl  =  G1G2F(31  ) 

GO  TO  100 
70  Gl  =  0. 

HI  a  0. 

GO  TO  100 
80  Gl  =  G1G2F (bl  ) 

90  HI  a  HPLMIF(B1.ZETA(J1),K) 

100  L  =  A 

ZETACa  ZETA(J2) 

X3  a  I.EIO 
Y3  =  l.ElO 
Z3  =  l.ElO 
DO  300  1=1.  L 

ANEW  =  (Y(J1  )-Y(J2)+X(J2)»ZETAC  -X ( J 1 ) *A1 ) / ( ZETAC  -Al) 
YNEW  =  Y( J1 ) +Al*(XNEw-X( J1 ) ) 

ZNEW  =  ZETA(  J1 )  +  ( (-1. )*»X*F1*(BETA(J2)-BETA(J1 ))  +  (-!. )*»X»G1» 
1  (  S(  J2  )-S(  J1  )  )-(-!  .  )»»X«Hl*(YNEii»-Y(  Jl)  )  /  Y(  J1  )  )/El 


IF 

(ABSF( (ZNEW-Z3)/ZNtW)-EPS) 

150. 

150. 

210 

150 

IF 

(ABSF( (XNEW-X3)/XNEW)-EPS) 

200. 

200. 

210 

200 

IF 

(ABSF( (YNEW-Y3)/YNEw)-EPS) 

500* 

500. 

210 

210 

X3 

=  ANEW 

Y3  =  YNEW 
Z3  =  ZNEW 

Z1=(ZETA( J1 )+ZNEW)/2. 
A1=APLMIF(B1 .Z1 *X) 
E1=E1E2F(Z1 ) 

IF  (HI)  260.  250.  260 
250  HI  =  C. 


40 


GO  TO  270 

260  =  )  *Y  (  J1 »  »2  .  /  <  Y  (  J1 ) +Y'NEW  ) 

270  2ETAC=<ZETA( J2)+ZNEW)/2. 

300  CONTINUE 

SENSE  LIGHT  4 
500  ZETA(J2)  =  ZNEW 
X{J2)  =  XNEW 
Y(J2)  =  YNEW 
RETURN 

1000  WRITE  OUTPUT  TAPE  6.  1.  M 
CALL  EXIT 

1  F0RMAT(38H1ERR0R  EXIT  FROM  FREE  SUBROUTINE.  M  =  16) 
END 
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Fig.  1  Illustration  for  the  method  of  computing 
and  storing  points  of  the  Mach  net. 


2  Illastration  for  the  method  of  computing 
new  point  along  a  shock. 


3  Illustration  for  method  of  reducing  net 
size  in  the  region  near  a  shock. 
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Fig.  4  Illastratlon  for  method  of  computing 

axially  symmetric  perfect  nozzles  using 
conical  source  flow  (Ref. [2]) 
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Fig.  5  Illustration  for  strip  method  of  computing 
Mach  net  of  the  tramsltion  region  of  an 
axially  symmetric  nozzle. 
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Fig.  6  Illustration  for  the  computation  of  the  wall 
streamline  through  the  transition  region  of 
an  axially  symmetric  nozzle. 
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ilg.  8  Illustration  for  the  calculation  of  the  wall 
streamline  through  the  Prandtl-Neyer  region 
of  a  two  dimensional  nozzle. 
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